      SUBROUTINE  COBSWF(CH,XI,M,N,MINT,INDX,NMAX,MAXJW,JW,D,ALAMD,     00010000
     1                   SIGMA,DSUM,R1,DR1,R2,DR2)                      00020000
CCC       THE OBLATE SPHEROIDAL RADIAL FUNCTIONS WITH COMPLEX ARGUMENTS.00030000
CCC       THE SIZE PARAMETER C IS COMPLEX.                              00040000
CCC       THE FIRST KIND FUNCTIONS AS A SERIES OF THE SPHERICAL BESSEL F00050000
CCC       INDX=1, EXTERNAL FIELD OR C1 ;  INDX=2, INTERNAL FIELD OR C2. 00060000
CCC       MINT=INDEX PARAMETER FOR THE CALCULATION OF BESSEL FUNCT.     00070000
CCC       THE SECOND KIND FUNCTIONS ARE NEVER CALCULATED FOR INDX = 2.  00080000
CCC                                                                     00090000
      PARAMETER  (JMX=150,MXPQ=190,MXSJ=150,MXSN=190)                   00100010
      COMPLEX*16 CI,C,C2,CH,CH2,XX,DX,W,T1,T2,T3,RATIO,SIGMA,PI,DSUM,   00110000
     1           ALAMD,AINT,E,F,G,S1,S2,S3,R1A,DR1A,RCONV1,R1,DR1,R2A,  00120000
     2           DR2A,RCONV2,R2N,DR2N,QSUM,DQSUM,QD,DN0,QNA,QNB,DQN,DP0,00130000
     3           PSUM,DPSUM,COEF2,R2Q,DR2Q,R2,DR2,HSUM,DHSUM,B,DWRAT,R3,00140000
     4           DR3,HCONV,R2H,DR2H,DR2NC,DR2QC,DR2HC                   00150000
      COMPLEX*16 D(JMX),DP(JMX),DN(80),DNRAT(80),DPRAT(JMX),PX(MXPQ),   00160010
     5           DPX(MXPQ),QX(MXPQ),DQX(MXPQ),QN(80),SJ(MXSJ,2),        00170010
     6           DSJ(MXSJ,2),SN(MXSN,2),DSN(MXSN,2),UP(70),DW(500),     00180003
     7           A(100),R2M(3),DR2M(3),AMX(JMX)                         00190006
      REAL*8  XI,XI2,XIA,XIB,ALEG,SJFACT,FACT,Q,FF,PAI2,RL,DRL,         00200012
     +        CRN,CRQ,CRH,FLOAT,SCLSN,SMLST,CMETHD,CR(4),               00210012
     +        RFAC1,RFAC2,RFAC3,RFAC4                                   00211012
      FLOAT(N)=DFLOAT(N)                                                00220000
      E(J)=C2*DFLOAT((M+M+J)*(M+M+J-1))/DFLOAT((2*(M+J)-1)*(2*(M+J)+1)) 00230000
      F(J)=DFLOAT((M+J)*(M+J+1))+C2*DFLOAT(2*(M+J)*(M+J+1)-2*M*M-1)/    00240000
     +     DFLOAT((2*(M+J)-1)*(2*(M+J)+3))-ALAMD                        00250003
      G(J)=C2*DFLOAT((J+1)*(J+2))/DFLOAT((2*(M+J)+1)*(2*(M+J)+3))       00260000
C                                                                       00270000
      PAI2=3.141592653589793D0/2.D0                                     00280003
      CI=(0.D0,1.D0)                                                    00290000
      C=(0.D0,-1.D0)*CH                                                 00300000
      C2= C*C                                                           00310003
      CH2=CH*CH                                                         00320000
      XI2=XI*XI+1.D0                                                    00330000
      XX=CI*XI                                                          00340000
      H=CDABS(CH)                                                       00350000
      HXX=H*XI*XI                                                       00360000
C                                                                       00370012
        IF(M.EQ.0)  THEN                                                00380003
      XIA=1.D0                                                          00390003
      XIB=0.D0                                                          00400003
        ELSE                                                            00410003
      XIA=DSQRT(XI2/(XI*XI))**M                                         00420000
      XIB=-DFLOAT(M)/(XI*XI2)                                           00430000
        END IF                                                          00440003
      L=N-M                                                             00450003
      PSCRT=(1.0E-12)*(10.0**(L/10.0))                                  00460003
      IF(MAXJW.GE.JMX-1) MAXJW= JMX-2                                   00470010
CC          CRITERION FOR THE METHOD OF COMPUTATION                     00480003
      IF(XI.GE.0.80) GO TO 5000                                         00490000
      ALP=0.00667*H+0.7334                                              00500000
      LM0=0.4*H-0.999                                                   00510000
      IF(XI.GE.0.5) LM0=0.6*H-1.999                                     00520000
      GO TO 5010                                                        00530000
 5000 LM0=H-0.999                                                       00540000
      IF(H.GT.25.0) LM0=1.9*H-23.499                                    00550000
      ALP=0.00667*H+0.73334                                             00560000
      IF(H.GT.28.0) ALP=1.3                                             00570000
 5010 MTHD=LM0-ALP*FLOAT(M)                                             00580000
      IF (MTHD.LT.0) MTHD=0                                             00590000
      LM1=MTHD-4                                                        00600000
      LM2=MTHD+4                                                        00610000
      LMTN=LM2+4                                                        00620000
      IF(H.GE.20.0.AND.M.GE.20) LMTN=LM2+FLOAT(M)*H/25.0                00630000
      CR(1)= 1.D+25                                                     00640000
      CR(2)= 1.D+10                                                     00650000
      CR(3)= 1.D+15                                                     00660000
      CR(4)= 1.D+30                                                     00670000
      R2=(0.D0,0.D0)                                                    00680000
      DR2=(0.D0,0.D0)                                                   00690000
CC        SPHERICAL BESSEL FUNCTIONS                                    00700003
      IF(M.EQ.MINT.AND.N.EQ.M) GO TO 4002                               00710000
      GO TO 4005                                                        00720000
 4002 NSJ=MAXJW+M                                                       00730000
       IF(NSJ.GE.MXSJ-1)  NSJ= MXSJ-2                                   00740011
      DX=CH*XI                                                          00750000
      W=(CDSIN(DX)-DX*CDCOS(DX))/(DX*DX)                                00760000
      T3=(0.0D0,0.0D0)                                                  00770000
      T2=(1.D-70,1.D-70)                                                00780000
      IF(DABS(DIMAG(DX)).LT.1.0D-10) T2=(1.D-70,0.D0)                   00790000
      NUP=10+DX/5.                                                      00800000
      NSWT=NSJ                                                          00810000
      NA=NSJ+NUP                                                        00820000
      DO 4000 NB=1,NA                                                   00830000
      NC=NA+1-NB                                                        00840000
      T1=DFLOAT(NC+NC+1)*T2/DX-T3                                       00850000
      IF(CDABS(T1)-1.D+70) 4001,4003,4003                               00860010
 4003 T1=T1*1.0D-70                                                     00870000
      T2=T2*1.0D-70                                                     00880000
      NSWT=NC                                                           00890000
 4001 T3=T2                                                             00900000
      T2=T1                                                             00910000
      IF(NC.GT.NSJ) GO TO 4000                                          00920000
      SJ(NC,INDX)=T1                                                    00930000
 4000 CONTINUE                                                          00940000
      RATIO=W/SJ(2,INDX)                                                00950000
      DO 4010 ND=1,NSJ                                                  00960000
      IF(ND.GT.NSWT) GO TO 4015                                         00970000
      SJ(ND,INDX)=RATIO*SJ(ND,INDX)                                     00980000
      GO TO 4010                                                        00990000
 4015 SJFACT=1.D-70                                                     01000000
      IF(CDABS(SJ(ND,INDX)).LT.1.0D-5) SJFACT=0.D0                      01010000
      SJ(ND,INDX)=RATIO*SJ(ND,INDX)*SJFACT                              01020000
 4010 CONTINUE                                                          01030000
      DSJ(1,INDX)=-SJ(2,INDX)                                           01040000
      DO 4020 NE=2,NSJ                                                  01050000
      NF=NE-1                                                           01060000
 4020 DSJ(NF+1,INDX)=SJ(NF,INDX)-DFLOAT(NF+1)*SJ(NF+1,INDX)/DX          01070000
CC        SPHERICAL NEUMANN FUNCTIONS                                   01080003
      IF(INDX.EQ.2) GO TO 4005                                          01090000
      IF(XI.LT.1.005) GO TO 4005                                        01100000
      IF(HXX.LT.15.0) GO TO 4005                                        01110000
      SCLSN=1.0E-20                                                     01120000
      IF(CDABS(DX).GE.15.0) SCLSN=1.0E-10                               01130000
      SN(1,INDX)=-CDCOS(DX)/DX                                          01140000
      SN(2,INDX)=-(CDCOS(DX)+DX*CDSIN(DX))/(DX*DX)                      01150000
      SN(1,INDX)=SN(1,INDX)*SCLSN                                       01160000
      SN(2,INDX)=SN(2,INDX)*SCLSN                                       01170000
      DSN(1,INDX)=-SN(2,INDX)                                           01180000
      MAXSN=MAXJW+M                                                     01190000
       IF(MAXSN.GE.MXSN-1) MAXSN= MXSN-2                                01200010
      DO 4030 NP=2,MAXSN                                                01210000
      NQ=NP-1                                                           01220000
      IF(NQ.LE.1) GO TO 4040                                            01230000
      SN(NQ+1,INDX)=DFLOAT(NQ+NQ-1)*SN(NQ,INDX)/DX-SN(NQ-1,INDX)        01240000
 4040 DSN(NQ+1,INDX)=SN(NQ,INDX)-DFLOAT(NQ+1)*SN(NQ+1,INDX)/DX          01250000
 4030 CONTINUE                                                          01260000
 4005 CONTINUE                                                          01270000
C         ASSOCIATED LEGNDRE FUNCTIONS                                  01280000
      IF(INDX.EQ.2) GO TO 4007                                          01290000
      IF(N.GT.M) GO TO 4007                                             01300000
      MAXPQ=MAXJW+M                                                     01310000
       IF(MAXPQ.GE.MXPQ-1)  MAXPQ= MXPQ-2                               01320010
      MLEG=0                                                            01330000
      ALEG=1.D0                                                         01340000
 4070 MLEG=MLEG+1                                                       01350000
      ALEG=ALEG*DFLOAT(2*MLEG-1)                                        01360000
      IF(MLEG.LT.M)  GO TO 4070                                         01370000
C                                                                       01380003
         CALL  CLEGQ(XX,M,MAXPQ,QX)                                     01390003
C                                                                       01400003
      PX(1)=ALEG*(CDSQRT(XX*XX-1.D0)**M)                                01410000
      PX(2)=DFLOAT(2*M+1)*XX*PX(1)                                      01420000
      DO 4050 I=1,MAXPQ                                                 01430000
      MI=M+I                                                            01440000
      IN=I-1                                                            01450000
      PX(I+2)=(DFLOAT(2*MI+1)*XX*PX(I+1)-DFLOAT(M+MI)*PX(I))/DFLOAT(I+1)01460000
      DPX(I)=(DFLOAT(I)*PX(I+1)-DFLOAT(M+I)*XX*PX(I))/(XX*XX-1.0)       01470000
      DPX(I)=CI*DPX(I)                                                  01480000
      IF(I.EQ.1) GO TO 4050                                             01490000
      DQX(I-1)=(DFLOAT(IN-M)*QX(I)-DFLOAT(IN)*XX*QX(I-1))/(XX*XX-1.0)   01500000
      DQX(I-1)=CI*DQX(I-1)                                              01510000
 4050 CONTINUE                                                          01520000
      IF(M.EQ.0) GO TO 4007                                             01530000
      S1=QX(2)                                                          01540000
      S2=QX(1)                                                          01550000
      DO 4080 KA=1,M                                                    01560000
      KK=KA-1                                                           01570000
      S3=(DFLOAT(1-2*KK)*XX*S2+DFLOAT(KK+M-1)*S1)/DFLOAT(M-KK)          01580000
      S1=S2                                                             01590000
      S2=S3                                                             01600000
      QN(KA)=S3                                                         01610000
 4080 CONTINUE                                                          01620000
 4007 CONTINUE                                                          01630000
CC        EIGENVALUES AND EXPANSION COEFFICIENTS - D                    01640000
      IF(N.GT.M) GO TO 4009                                             01650000
         CALL  CMATR(M,CH,NMAX,AMX,-1)                                  01660003
 4009 MTRX=L+1                                                          01670000
      AINT=AMX(MTRX)                                                    01680000
C                                                                       01690011
         CALL  CDCOEF(CH,M,N,-1,AINT,JW,D,ALAMD,PI,SIGMA,DSUM,FCTR)     01700000
C                                                                       01710003
      IF(CDABS(DSUM).LT.1.0D-50) GO TO 3                                01720000
CC            RADIAL FUNCTION AND ITS DERIVATIVE - RJ(M,N,HX),DRJ       01730003
CC            THE RADIAL FUNCTIONS OF FIRST KIND                        01740003
      IFIN=JW                                                           01750000
      IF(MOD(L,2).EQ.0) IFIN=JW+1                                       01760000
      R1A=(0.D0,0.D0)                                                   01770000
      DR1A=(0.D0,0.D0)                                                  01780000
      DO 1060 JR=1,IFIN,2                                               01790000
      J=JR-1                                                            01800000
      IF(MOD(L,2).EQ.1) J=JR                                            01810000
      MR=M+J                                                            01820000
      IF(MR.GT.NSJ) GO TO 1060                                          01830000
      IR=IABS(J-L)                                                      01840000
      Q=-1.D0                                                           01850000
      IF(MOD(IR,4).EQ.0) Q=1.D0                                         01860000
      RFAC1= 1.D0                                                       01870011
      RFAC2= 1.D0                                                       01880011
        IF(M.GE.1)  THEN                                                01890011
      DO 15 IP=1, M                                                     01900011
      RFAC1= RFAC1*DFLOAT(J+IP)                                         01910011
   15 RFAC2= RFAC2*DFLOAT(J+IP+M)                                       01920011
        END IF                                                          01930011
      RCONV1= D(JR)*RFAC1*DSJ(MR+1,INDX)*RFAC2/FCTR                     01940013
      R1A=R1A+Q*D(JR)*RFAC1*SJ(MR+1,INDX)*RFAC2/FCTR                    01950013
      DR1A=DR1A+Q*RCONV1                                                01960000
      IF(CDABS(RCONV1/DR1A).LE.1.0D-15) GO TO 5                         01970000
 1060 CONTINUE                                                          01980000
    5 R1=XIA*R1A/PI                                                     01990013
      DR1=XIB*R1+CH*XIA*DR1A/PI                                         02000013
      IF(INDX.EQ.2) GO TO 4                                             02010000
CC         RADIAL FUNCTION OF SECOND KIND IN TERMS OF SPHERICAL NEUMANN 02020003
      IF(XI.LT.1.005) GO TO 7                                           02030000
      IF(HXX.LT.15.0) GO TO 7                                           02040000
      IF(L.GT.LMTN) GO TO 7                                             02050000
      R2A=(0.D0,0.D0)                                                   02060000
      DR2A=(0.D0,0.D0)                                                  02070000
      DO 1070 JR=1,IFIN,2                                               02080000
      J=JR-1                                                            02090000
      IF(MOD(L,2).EQ.1) J=JR                                            02100000
      MR=M+J                                                            02110000
      IF(MR.GE.MAXSN) GO TO 1070                                        02120000
      IR=IABS(J-L)                                                      02130000
      Q=-1.                                                             02140000
      IF(MOD(IR,4).EQ.0) Q=1.                                           02150000
      RFAC3= 1.D0                                                       02160013
      RFAC4= 1.D0                                                       02170011
        IF(M.GE.1)  THEN                                                02180011
      DO 17 IP=1, M                                                     02190011
      RFAC3=RFAC3*DFLOAT(J+IP)                                          02200011
   17 RFAC4=RFAC4*DFLOAT(J+IP+M)                                        02210011
        END IF                                                          02220011
      RCONV2= D(JR)*RFAC3*SN(MR+1,INDX)*RFAC4/FCTR               9      02230013
      R2A=R2A+Q*RCONV2                                                  02240000
      DR2A=DR2A+Q*D(JR)*RFAC3*DSN(MR+1,INDX)*RFAC4/FCTR                 02250013
      IF(CDABS(RCONV2/R2A).LE.1.0D-08) GO TO 8                          02260000
 1070 CONTINUE                                                          02270000
    8 R2N=XIA*R2A/PI                                                    02280013
      DR2N=XIB*R2N+CH*XIA*DR2A/PI                                       02290013
      R2N=R2N/SCLSN                                                     02300000
      DR2N=DR2N/SCLSN                                                   02310000
      DR2NC=(1.D0/(CH*XI2)+R2N*DR1)/R1                                  02320000
      CRN=CDABS((DR2NC-DR2N)/DR2NC)                                     02330000
      R2M(1)= R2N                                                       02340000
      DR2M(1)=DR2N                                                      02350000
      CR(1)= CRN                                                        02360000
    7 CONTINUE                                                          02370000
CCC       SPHEROIDAL RADIAL FUNCTIONS OF THE SECOND KIND                02380000
CCC       IN TERMS OF THE ASSOCIATED LEGENDRE FUNCTIONS                 02390000
      IF(L.LT.LM1) GO TO 1050                                           02400000
      IF(MOD(L,2).EQ.1) GO TO 3000                                      02410000
      JFIN=2*M                                                          02420000
      JINT=2                                                            02430000
      GO TO 3010                                                        02440000
 3000 JFIN=2*M-1                                                        02450000
      JINT=1                                                            02460000
 3010 QSUM=(0.D0,0.D0)                                                  02470000
      DQSUM=(0.D0,0.D0)                                                 02480000
      DO 10 I=1,IFIN,2                                                  02490000
      IA=I-1                                                            02500000
      IF(MOD(L,2).EQ.1)  IA=I                                           02510000
      MA=M+IA                                                           02520000
      IF(MA.GE.MAXPQ) GO TO 11                                          02530000
      QD=D(I)*DQX(MA+1)                                                 02540000
      QSUM=QSUM+D(I)*QX(MA+1)                                           02550000
      DQSUM=DQSUM+QD                                                    02560000
      IF(CDABS(QD/DQSUM).LT.1.0D-13) GO TO 11                           02570000
   10 CONTINUE                                                          02580000
   11 IF(M.EQ.0) GO TO 3030                                             02590000
      DNRAT(JFIN+2)=(0.D0,0.D0)                                         02600000
      DO 20 J=JINT,JFIN,2                                               02610000
      JA=JINT+JFIN-J                                                    02620000
   20 DNRAT(JA)=-E(-JA+2)/(F(-JA)+G(-JA-2)*DNRAT(JA+2))                 02630000
      DN0=D(1)                                                          02640000
      DO 30 JB=JINT,JFIN,2                                              02650000
      DN(JB)=DN0*DNRAT(JB)                                              02660000
      MJB=M-JB                                                          02670000
      MJC=MJB+1                                                         02680000
      IF(MJB.LT.0) GO TO 2050                                           02690000
      QNA=QX(MJB+1)                                                     02700000
      QNB=QX(MJC+1)                                                     02710000
      GO TO 2060                                                        02720000
 2050 MJD=-MJB                                                          02730000
      IF(MJD.EQ.1) GO TO 2070                                           02740000
      QNA=QN(MJD)                                                       02750000
      QNB=QN(MJD-1)                                                     02760000
      GO TO 2060                                                        02770000
 2070 QNA=QN(1)                                                         02780000
      QNB=QX(1)                                                         02790000
 2060 DQN=(FLOAT(MJB-M+1)*QNB-FLOAT(MJB+1)*XX*QNA)/(XX*XX-1.0)          02800000
      DQN=CI*DQN                                                        02810000
      QSUM=QSUM+DN(JB)*QNA                                              02820000
      DQSUM=DQSUM+DN(JB)*DQN                                            02830000
      DN0=DN(JB)                                                        02840000
   30 CONTINUE                                                          02850000
      IF(MOD(L,2).EQ.1) GO TO 3040                                      02860000
      COEF2=FF(M,N)*FACTMM(M)*DN(JFIN)*C*PI/(DFLOAT(M+M-1)*C**M)*       02870011
     +      FCTR                                                        02880011
      GO TO 3050                                                        02890000
 3040 COEF2=-FF(M,N)*FACTMM(M)*DN(JFIN)*C2*PI/(DFLOAT((M+M-3)*(M+M-1))* 02900011
     +      C**M)*FCTR                                                  02910011
 3050 CONTINUE                                                          02920000
      GO TO 3060                                                        02930000
 3030 IF(MOD(L,2).EQ.1) GO TO 3070                                      02940000
      COEF2=-FF(M,N)*D(1)*C*PI                                          02950000
      GO TO 3060                                                        02960000
 3070 COEF2=-FF(M,N)*D(1)*C2*PI/3.0D0                                   02970000
 3060 CONTINUE                                                          02980000
      IH=H                                                              02990000
      IF(MOD(L,2).EQ.1) GO TO 3080                                      03000000
      KINT=M+M+2                                                        03010000
      KFIN=30+2*(M+IH)                                                  03020000
      GO TO 3090                                                        03030000
 3080 KINT=M+M+1                                                        03040000
      KFIN=31+2*(M+IH)                                                  03050000
 3090 IF(KFIN.GE.138) KFIN=138+MOD(L,2)                                 03060010
 2040 DPRAT(KFIN+2)=(0.D0,0.D0)                                         03070000
      DP0=D(1)                                                          03080000
      IF(M.NE.0)  DP0=DN(JFIN)                                          03090000
      DO 40 KA=KINT,KFIN,2                                              03100000
      KB=KINT+KFIN-KA                                                   03110000
      IF(KB.EQ.KINT) GO TO 2010                                         03120000
      DPRAT(KB)=-E(-KB+2)/(F(-KB)+G(-KB-2)*DPRAT(KB+2))                 03130000
      GO TO 40                                                          03140000
 2010 IF(MOD(L,2).EQ.1) GO TO 2020                                      03150000
      DPRAT(KB)=C2/FLOAT((M+M+1)*(M+M-1))/(F(-KB)+G(-KB-2)*DPRAT(KB+2)) 03160000
      GO TO 40                                                          03170000
 2020 DPRAT(KB)=-C2/FLOAT((M+M-1)*(M+M-3))/(F(-KB)+G(-KB-2)*DPRAT(KB+2))03180000
   40 CONTINUE                                                          03190000
      PSUM=(0.D0,0.D0)                                                  03200000
      DPSUM=(0.D0,0.D0)                                                 03210000
      DO 50 KC=KINT,KFIN,2                                              03220000
      DP(KC)=DP0*DPRAT(KC)                                              03230000
      KD=KC-M-1                                                         03240000
      KE=KD+1-M                                                         03250000
      IF(KE.GT.MAXPQ) GO TO 50                                          03260000
      PSUM=PSUM+DP(KC)*PX(KE)                                           03270000
      DPSUM=DPSUM+DP(KC)*DPX(KE)                                        03280000
      DP0=DP(KC)                                                        03290000
   50 CONTINUE                                                          03300000
      KFINM=KFIN-M-M                                                    03310000
      IF(KFINM.GE.MAXPQ) GO TO 2030                                     03320000
        IF(KFIN.GE.137) GO TO 2030                                      03330010
      CRIT=CDABS(DP(KFIN)*PX(KFINM)/PSUM)                               03340000
      KFIN=KFIN+4                                                       03350000
      IF(CRIT.GT.PSCRT) GO TO 2040                                      03360000
 2030 KFIN=KFIN-4                                                       03370000
      R2Q=(QSUM+PSUM)/COEF2                                             03380000
      DR2Q=(DQSUM+DPSUM)/COEF2                                          03390000
      DR2QC=(1.D0/(CH*XI2)+R2Q*DR1)/R1                                  03400000
      CRQ=CDABS((DR2QC-DR2Q)/DR2QC)                                     03410000
      R2M(2)= R2Q                                                       03420000
      DR2M(2)=DR2Q                                                      03430000
      CR(2)= CRQ                                                        03440000
 1050 CONTINUE                                                          03450000
C            THE RADIAL FUNCTION OF THIRD KIND BY HANISH,S METHOD       03460000
      IF(H.LT.7.0) GO TO 1080                                           03470000
      IF (L.GT.LMTN) GO TO 1080                                         03480000
      LMT=3*(H+CDABS(ALAMD)/H)                                          03490000
      NUP=LMT+20                                                        03500000
      ISTP=MAXPQ-2                                                      03510000
      IF(NUP.LE.ISTP) NUP=ISTP+10                                       03520000
      IF(NUP.GE.499) NUP=499                                            03530000
      DO 55 I=1,ISTP                                                    03540000
   55 A(I)=(0.D0,0.D0)                                                  03550000
      MP1=M+L+1                                                         03560000
      UP(1)=-3.0*(ALAMD+CH2)/(2.0*CH*DFLOAT(M+1))                       03570000
      DO 60 IR=1,MP1                                                    03580000
      I=IR-M                                                            03590000
      MI=M+I                                                            03600000
      UP(IR+1)=DFLOAT(2*MI+3)*(DFLOAT(MI*(MI+1))-ALAMD-CH2) /(2.0*CH*   03610000
     +  DFLOAT((MI+1)*(MI+M+1)))+DFLOAT(2*MI+3)*DFLOAT(I*MI)/(UP(IR)*   03620007
     +  DFLOAT(2*MI-1)*DFLOAT((MI+1)*(MI+M+1)))                         03630007
   60 CONTINUE                                                          03640000
      A(1)=UP(1)                                                        03650000
      DO 62 I=1,MP1                                                     03660000
   62 A(I+1)=A(I)*UP(I+1)                                               03670000
      DW(NUP+1)=(CH+ALAMD/CH)/DFLOAT(NUP-M)                             03680000
      DO 64 JR=MP1,NUP                                                  03690000
      J=(NUP+MP1)-JR-M                                                  03700000
      MJ=M+J                                                            03710000
      IF(J.LE.0) GO TO 64                                               03720000
      DW(MJ)=2.0*CH*DFLOAT(J*MJ)/DFLOAT(2*MJ-1)/(2.0*CH*DFLOAT((MJ+1)*  03730000
     + (MJ+M+1))*DW(MJ+1)/DFLOAT(2*MJ+3)-(DFLOAT(MJ*(MJ+1))-ALAMD-CH2)) 03740008
   64 CONTINUE                                                          03750000
      DWRAT=UP(MP1)/DW(MP1)                                             03760000
      ISTP1=ISTP+1                                                      03770000
      DO 66 I=MP1,ISTP1                                                 03780000
   66 DW(I)=DWRAT*DW(I)                                                 03790000
      DO 68 K=MP1,ISTP                                                  03800000
      IF(CDABS(A(K)).LT.1.0D-50) GO TO 69                               03810000
      IF(CDABS(A(K)).GT.1.D+70) GO TO 69                                03820000
   68 A(K+1)=A(K)*DW(K+1)                                               03830000
   69 CONTINUE                                                          03840000
      B=CI**(M+M+1)/(CH*FACT(M))*CDEXP(CI*(CH*XI-DFLOAT(N+1)*PAI2))     03850000
      HSUM=QX(1)                                                        03860000
      DHSUM=DQX(1)                                                      03870000
      DO 70 I=1,ISTP                                                    03880000
      HCONV=A(I)*DQX(I+1)                                               03890000
      HSUM=HSUM+A(I)*QX(I+1)                                            03900000
      DHSUM=DHSUM+HCONV                                                 03910000
      IF(CDABS(HCONV/DHSUM).LE.1.D-18) GO TO 72                         03920000
   70 CONTINUE                                                          03930000
   72 R3=B*HSUM                                                         03940000
      DR3=CI*CH*R3+B*DHSUM                                              03950000
      R2H=(R3-R1)/CI                                                    03960000
      DR2H=(DR3-DR1)/CI                                                 03970000
      IF(DABS(DIMAG(CH)).GT.1.D-05) GO TO 1                             03980000
      RL=R2H                                                            03990000
      DRL=DR2H                                                          04000000
      R2H=(1.D0,0.D0)*RL                                                04010000
      DR2H=(1.D0,0.D0)*DRL                                              04020000
    1 DR2HC=(1.D0/(CH*XI2)+R2H*DR1)/R1                                  04030000
      CRH=CDABS((DR2HC-DR2H)/DR2HC)                                     04040000
      R2M(3)= R2H                                                       04050000
      DR2M(3)=DR2H                                                      04060000
      CR(3)= CRH                                                        04070000
 1080 CONTINUE                                                          04080000
C         CHOOSSING THE BEST VALUES.                                    04090000
      SMLST=DMIN1(CR(1),CR(2),CR(3))                                    04100008
      DO 7000 II=1,3                                                    04110000
      I=II                                                              04120000
      IF(CR(II).LT.1.D-13) GO TO 7010                                   04130000
      CMETHD=DABS((SMLST-CR(II))/CR(II))                                04140000
      IF(CMETHD.LT.1.D-7) GO TO 7010                                    04150000
 7000 CONTINUE                                                          04160000
      R2= R2M(2)                                                        04170000
      DR2=DR2M(2)                                                       04180000
      IF(L.GT.MTHD) GO TO 7020                                          04190000
      R2= R2M(3)                                                        04200000
      DR2=DR2M(3)                                                       04210000
      GO TO 7020                                                        04220000
 7010 R2= R2M(I)                                                        04230000
      DR2=DR2M(I)                                                       04240000
 7020 CONTINUE                                                          04250000
      GO TO 4                                                           04260000
    3 WRITE(6,100)                                                      04270000
  100 FORMAT(1H0,20X,'*  CPRSWF  -  RETURN WITHOUT CALCULATION  *  COEFF04280000
     1ICIENTS D ARE ALL ZERO' )                                         04290000
    4 RETURN                                                            04300000
      END                                                               04310000
